Задание 1¶
Сотрудник службы оценки качества продукции заметил, что число бракованных листов стали в партии для стали марки A больше, чем для стали марки B. Также он обратил внимание, что при средней скорости прокатки более 4 м/с число бракованных листов больше. В соответствии с этими наблюдениями, предлагается снизить скорость прокатки и ввести дополнительные меры контроля качества для стали марки A.
Обоснуйте, что:
1) Более 3 бракованных листов на партию выходит значимо чаще для стали марки А, чем для стали марки B.
2) При скоростях прокатки более 4 м/с свыше 3 бракованных листов стали на партию выходит значимо чаще, чем при меньших скоростях прокатки.
Данные, необходимые для анализа, содержатся в файле «Статистика за 2018 год». Количество листов стали в каждой партии предполагается одинаковым.
import pandas as pd
import numpy as np
from statsmodels.stats.weightstats import _tconfint_generic
# считываем данные
data = pd.read_csv('Статистика за 2018 год.csv',encoding='windows-1251')
data
Решение¶
Для начала посмотрим на заполненность нашего датасета и распределение данных.
# количество пропусков в каждой категории
data.isnull().sum()
# pip install plotly
import plotly
import plotly.express as px
# распределение Числа бракованных листов в зависимости от Скорости прокатки
fig1 = px.histogram(data, y="Число бракованных листов", x = "Скорость прокатки", title = 'Распределение числа бракованных листов в зависимости от скорости прокатки')
fig1.show()
fig2 = px.histogram(data, x="Число бракованных листов",color_discrete_sequence=['indianred'], title = 'Распределение частоты количества числа бракованных листов в партии')
fig2.show()
Воспользуемся Критерием Шапиро-Уилка для оценки нормальности распределения числа бракованных листов:
$H_0\colon$ число бракованных листов распределено нормально
$H_1\colon$ число бракованных листов распределено не нормально.
import scipy
from scipy import stats
print("Тест Шапиро-Уилка W-statistic: %f, p-value: %f" % stats.shapiro(data['Число бракованных листов']))
Значение p-value < 0.05, следовательно, гипотеза H0 отвергается.
Вывод:¶
- Датасет заполнен корректно. Пропусков в данных нет.
- Критерий Шапиро-Уилка показал, что данные по числу бракованных листов не распределены нормально.
Вернемся к гипотезам, которые предлагаются сотрудник службы оценки качества.
- Более 3 бракованных листов на партию выходит значимо чаще для стали марки А, чем для стали марки B.
# определим для каждой марки стали датасет с числом бракованных листов более 3 на партию
# марка стали А
data_A = data[(data['Марка стали']=='A') & (data['Число бракованных листов']>3)]['Число бракованных листов']
# марка стали B
data_B = data[(data['Марка стали']=='B') & (data['Число бракованных листов']>3)]['Число бракованных листов']
Посмотрим на распределения по каждой марки стали.
import plotly.express as px
# распределение Числа бракованных листов в зависимости от Скорости прокатки
fig1 = px.histogram(data_A, title = 'Марка стали "А"')
fig1.show()
fig2 = px.histogram(data_B,color_discrete_sequence=['indianred'], title = 'Марка стали "В"')
fig2.show()
Рассчитаем доверительные интервалы для оценки среднего. Поскольку у нас неизвестна дисперсия распределения, то используем выборочные дисперсии, и построим доверительные интервалы вида: $$\bar{X}_n \pm t_{1-\frac{\alpha}{2}} \frac{S}{\sqrt{n}}$$
# функция для расчета доверительных интервалов для оценки среднего
# _tconfint_generic подаются значения(выборочного среднего(Xn),выборочной дисперсии(S),параметр alpha,выбор альтернати)
def confidence_interval(df):
# 95% доверительный интервал
interval = _tconfint_generic(df.mean(),
df.std(ddof=1)/np.sqrt(df.shape[0]),
df.shape[0]-1, 0.05, 'two-sided')
return {'95% доверительный интервал': interval}
# значения для марки стали "А" и "B"
print('Значения для марки стали "А":', confidence_interval(data_A))
print('Значения для марки стали "В":', confidence_interval(data_B))
Доверительные интервалы пересекаются сделать вывод о том, что гипотеза верна сложно.
Ранговый критерий Манна-Уитни¶
Воспользуемся ранговым критерием Манна-Уитни для оценки различий между двумя выборками.
H0: Вероятность того, что значение из первой выборки будет больше, чем значение из второй выборки равна вероятности того, что значение из второй выборки будет больше, чем значение из первой выборки.
H1: это не так
stats.mannwhitneyu(data_A, data_B)
Поскольку значение p-value > 0.05 , то можно сделать вывод, что гипотеза H0 не отвергается.
Вывод:¶
По доверительному интервалу ожно заметить, что по марке стали "А" более 3 бракованных листов на партию встречается чаще, чем по марке стали "B".
Однако критерий Манна-Уитни показывает, что вероятность того, что значение из первой выборки будет больше, чем значение из второй выборки равна вероятности того, что значение из второй выборки будет больше, чем значение из первой выборки, то есть гипотезу о том, что "Более 3 бракованных листов на партию выходит значимо чаще для стали марки А, чем для стали марки B" нельзя считать правдивой.
- При скоростях прокатки более 4 м/с свыше 3 бракованных листов стали на партию выходит значимо чаще, чем при меньших скоростях прокатки.
# определим датасеты в зависимости от Скорости прокатки и количества бракованных листов
# скорость прокатки более 4 м/с
speed_more = data[(data['Скорость прокатки']>4) & (data['Число бракованных листов']>3)]['Число бракованных листов']
# скорость прокатки менее 4 м/с
speed_less = data[(data['Скорость прокатки']<=4) & (data['Число бракованных листов']>3)]['Число бракованных листов']
import plotly.express as px
# распределение Числа бракованных листов в зависимости от Скорости прокатки
fig1 = px.histogram(speed_more, title = 'Значения для скорости прокатки более 4 м/с')
fig1.show()
fig2 = px.histogram(speed_less,color_discrete_sequence=['indianred'], title = 'Значения для скорости прокатки менее 4 м/с')
fig2.show()
# значения для марки стали "А" и "B"
print('Значения для скорости прокатки более 4 м/с:', confidence_interval(speed_more))
print('Значения для скорости прокатки менее 4 м/с:', confidence_interval(speed_less))
Ранговый критерий Манна-Уитни¶
Воспользуемся ранговым критерием Манна-Уитни для оценки различий между двумя выборками.
H0: Вероятность того, что значение из первой выборки будет больше, чем значение из второй выборки равна вероятности того, что значение из второй выборки будет больше, чем значение из первой выборки.
H1: это не так
stats.mannwhitneyu(speed_more, speed_less)
Поскольку значение p-value > 0.05 , то можно сделать вывод, что гипотеза H0 не отвергается.
Вывод:¶
Доверительный интервал для значений скорости прокатки более 4 м/с находится внутри доверительного интервала для значений скорости прокатки менее 4 м/с.
Это означает, что для данного датасета при 95% довертительном интервале гипотеза о том, что "при скоростях прокатки более 4 м/с свыше 3 бракованных листов стали на партию выходит значимо чаще, чем при меньших скоростях прокатки", - не является верной.
Критерий Манна-Уитни для оценки различий между двумя выборками также подтверждает, что фактор скорости не влияет на увеличения количества бракованных листов свыше 3 на партию.
Задание 2¶
Вам необходимо построить модель, которая на основании данных, поступающих каждую минуту, определяют качество продукции, производимое на обжиговой машине.
Обжиговая машина представляет собой агрегат, состоящий из 5 одинаковых по размеру камер, в каждой камере установлено по 3 датчика температур. Кроме этого, для данной задачи Вы собрали данные о высоте слоя сырья и его влажности. Высота слоя и влажность измеряются при входе сырья в машину. Сырье проходит через обжиговую машину за час.
Данные с показателями работы обжиговой машины содержатся в файле X_data.csv:
Качество продукции измеряется в лаборатории по пробам, которые забираются каждый час, данные по известным анализам содержатся в файле Y_train.csv. В файле указано время забора пробы, проба забирается на выходе из обжиговой машины.
Вы договорились с заказчиком, что оценкой модели будет являться показатель MAE, для оценки модели необходимо сгенерировать предсказания за период, указанный в файле Y_submit.csv (5808 предиктов).
Все необходимые данные лежат по ссылке: https://www.dropbox.com/s/lo2w549fv8lo3oc/Test2.7z
Решение¶
# считаем данные
X_data = pd.read_csv('X_data.csv', sep=';') # данные с показателями работы обжиговой машины
Y_train = pd.read_csv('Y_train.csv',sep=';', header=None)
Y_submit = pd.read_csv('Y_submit.csv',sep=';', header=None)
# данные с показателями работы обжиговой машины
X_data.head(5)
# данные по известным анализам
#В файле указано время забора пробы, проба забирается на выходе из обжиговой машины.
Y_train.head(5)
# данные по которым необходимо сгенерировать предсказания за период
Y_submit.head(5)
1. Краткий анализ X_data.¶
# размер датасета
X_data.shape
X_data.info(null_counts=True)
Вывод:¶
В датасете нет пропущеных значений.
2. Анализ поставленной задачи.¶
Из условия задачи известно, что "Высота слоя (H_data) и влажность (AH_data) измеряются при входе сырья в машину. Сырье проходит через обжиговую машину за час".
Следовательно, можно сделать следующий вывод:
Сырье поступает в обжиговую машину во время, которое указано в первом столбце ('Unnamed: 0') X_data. Чтобы получить информацию по качеству продукции, необходимо к X_data['Unnamed: 0'] добавить 1 час и сопоставить эти данные с Y_train.
from datetime import datetime, timedelta
# добавим к времени входа + 1 час
X_data['Unnamed: 0'] = pd.to_datetime(X_data['Unnamed: 0']) + timedelta(hours=1)
X_data[['Unnamed: 0']]
# сопоставим значения из Y_train с X_data
data_for_modeling = X_data.loc[X_data['Unnamed: 0'].isin(Y_train[0])]
print('Размерность получившегося датасета: ', data_for_modeling.shape)
print('Размерность Y_train: ', Y_train.shape)
data_for_modeling
Посмотрим на распределение признаков в data_for_modeling
for i in list(data_for_modeling.loc[:, data_for_modeling.columns != 'Unnamed: 0'].columns):
fig1 = px.histogram(data_for_modeling[i], title = i)
fig1.show()
Вывод:¶
Сформировали датасет, который будет использоваться для моделирования, путем добавления времени обжига (1 час) ко времени входа сырья в машина и сопоставлением этого датафрейма с Y_train.
Получился датасет размером (29184, 18). В целом, признаки в датасете не имеют каких-то аномалий и значительных выбросов. Можно приступать к моделированию.
3. Разделение данных на train и test.¶
from sklearn.model_selection import cross_val_score, train_test_split
# убереме из датасет столбец с информцией о времени
X_train, X_test, y_train, y_test = train_test_split(data_for_modeling.loc[:, data_for_modeling.columns != 'Unnamed: 0'],
Y_train[1], shuffle=True,
test_size=0.2, random_state=3)
# размерность данных для обучения и тестирования
print('Данные для обучения: ', X_train.shape,y_train.shape)
print('Данные для теста: ',X_test.shape, y_test.shape)
Вывод:¶
Сформировали датасеты для обучения и теста в соотношении 80/20.
4. Построение модели.¶
Построим несколько моделей, чтобы сравнить какой алгоритм покажет себя лучше. Остановимся на моделях:
- Градиентоный бустинг с подбором параметров.
- Случайный лес с подборо параметров.
- Xgboost из коробки.
Подберем параметры, обучим алгоритм и оценим метрику MAE для каждого алгоритма.
4.1. Градиентный бустинг.¶
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import GridSearchCV
import warnings
warnings.filterwarnings('ignore')
Сделаем небольшой подбор гиперпараметров при помощи GridSearch
# список параметров для бустинга
boosting_params = {
'n_estimators': list(range(100,500,100)),
'min_samples_leaf': list(range(2,6,2)),
'max_depth': list(range(5,10,2)),
'max_features': list(range(5,10,2)),
}
boosting_gs = GridSearchCV(GradientBoostingRegressor(random_state=3),
boosting_params,
scoring='neg_mean_absolute_error',
n_jobs=-1,
verbose=True)
# обучаем GridSearchCV
boosting_gs.fit(X_train, y_train)
boosting_gs.best_params_
print('Параметры для алгоритма:')
print('max_depth: 9, max_features: 9, min_samples_leaf: 4, n_estimators: 200')
#укажем параметры в алгоритме
model_gb = GradientBoostingRegressor(max_depth=9, max_features= 9,
min_samples_leaf=4, n_estimators=200,
random_state=3)
# обучаем модель
model_gb.fit(X_train, y_train)
# делаем предсказания
predictions_model_gb = model_gb.predict(X_test)
from sklearn.metrics import mean_absolute_error
print('Значение метрики МАE на отложенной выборке для алгоритма градиетного бустинга:', mean_absolute_error(y_test, predictions_model_gb))
4.2. Случайные лес¶
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
import warnings
warnings.filterwarnings('ignore')
# список параметров для случайного леса
forest_params = {
'n_estimators': list(range(1500,2000, 100)),
'max_depth': list(range(22, 28, 2)),
'max_features': list(range(8,16,2)),
}
forest_gs = GridSearchCV(RandomForestRegressor(n_jobs=-1, random_state=3),
forest_params,
scoring='neg_mean_absolute_error',
n_jobs=-1,
verbose=True)
forest_gs.fit(X_train, y_train)
forest_gs.best_params_
print('Параметры для алгоритма:')
print('max_depth: 26, max_features: 12, n_estimators: 1900')
model_rf = RandomForestRegressor(max_depth=26, max_features=12, n_estimators=1900,
random_state=3)
# обучаем модель
model_rf.fit(X_train, y_train)
# делаем предсказания
predictions_model_rf = model_rf.predict(X_test)
from sklearn.metrics import mean_absolute_error
print('Значение метрики МАE на отложенной выборке для алгоритма случайный лес:', mean_absolute_error(y_test, predictions_model_rf))
4.3. Xgboost¶
import xgboost as xgb
model_xgb = xgb.XGBRegressor(random_state=3)
# обучаем модель
model_xgb.fit(X_train, y_train)
# делаем предсказания
predictions_model_xgb = model_xgb.predict(X_test)
from sklearn.metrics import mean_absolute_error
print('Значение метрики МАE на отложенной выборке для алгоритма xgboost:', mean_absolute_error(y_test, predictions_model_xgb))
Вывод:¶
Модель градиентного бустинга с подбором параметров показала себя лучше, поэтому остановимся на ней.
5. Вычисление качества продукции для Y_submit.¶
Y_submit.head(5)
# предсказания за период, указанный в файле Y_submit
X_submit = X_data.loc[X_data['Unnamed: 0'].isin(Y_submit[0])]
X_submit.shape, Y_submit.shape
X_submit
# генерируем предсказания в Y_pred
Y_pred = pd.DataFrame(model_gb.predict(X_submit.loc[:, X_submit.columns != 'Unnamed: 0']))
Y_pred
# выгрузка в csv
Y_pred.to_csv('Y_pred.csv')
Вывод:¶
Можно также постараться улучшить модель путем подбора более сложных алгоритмов (CatBoost, LightGBM, Нейронные сети) и генерацией новых фичей на основе исходных признаков (квадраты признаков, квадратные корни из признаков, логарфимы и т.д.)